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Abstract. We present observations of 15 massive cores in three different CS transitions from the FCRAO 14m and the KOSMA 
3m telescope. We derive physical parameters of these cores using different approaches to the line radiative transfer problem. 
The local radiative transfer approximations fail to provide reliable values except for the column densities. A self-consistent 
explanation of the observed line profiles is only possible when taking density gradients and an internal turbulent structure of 
the cores into account. The observational data can be fitted by a spherically symmetric radiative transfer model including such 
gradients and a turbulent clumping. We find that the observed massive cores are approximately virialised with a clumpy density 
profile that decays with a radial exponent of about —1.6 down to a relatively sharp outer boundary. 

We show that a careful analysis of spatially unresolved multi-line observations using a physical radiative transfer model can 
provide values for physical parameters that could be obtained otherwise only by direct observations with much higher spatial 
resolution. This applies to all quantities directly affecting the line excitation, like the mass and size of dense cores. Information 
on the exact location or number of clumps, of course, always has to rely on high-resolution observations e.g. from interferom- 
eters. 
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1. Introduction 

Whereas the average density in molecular clouds falls between 
about 50 and 1000 cm""^, observations of high-dipole-moment 
molecules like CS or NH3 reveal cores with densities up to 10® 
cm^"^. Massive cores are typically somewhat warmer than their 
embedding molecular cloud and show sizes of about 1 — 3 pc 
and masses between some ten and some thousand solar masses. 
They appear as relatively bright objects in molecular line maps 
so that they are favoured objects from the viewpoint of the ob- 
servations. The observed line widths of 2 — 15 km/s are consid- 
erably larger than thermal, indicating turbulent motions, pos- 
sibly induced by outflows from star-forming activity. Many of 
them are associated with young OB or T Tauri stars suggesting 
that massive cores are sites of massive and multiple star forma- 



tion (c.f. Myers 1999). For a better understanding of star for- 



mation we have to know the physical parameters within these 
cores, i.e. the geometrical structure, the density, velocity, and 
temperature distribution. 

Dense cores are best traced by molecules like CS, HC3N, 
NH3, and H2CO characterised by high critical densities even 
for low transitions. Throughout this paper we will use the term 
"core" for the dense inner part of a cloud which is directly visi- 
ble in CS, whereas the rest of the cloud may contribute to the ra- 
diative excitation but is mostly visible in CO rather than in CS. 
We call smaller substructures within the cores "clumps". Meta- 
stable transitions of NH3 provide reasonable estimates for the 



kinetic core temperature and the simultaneous observation of 
several lines from the same isotope provide combined informa- 
tion both on the density and the temperature structure. Here, 
we use observations of the 2-1, 5^, and 7-6 transitions of CS 
and C^^'S obtained with the FCRAO 14m and the KOSMA 3m 
telescopes to study 15 massive cores and complement our ob- 
servations with CS data from the literature to derive the physi- 
cal core parameters. 

Plume et al. (1997) observed about 150 massive cores in 



CS and C''"'S using the IRAM 30m telescope. With an escape 
probability approximation they derive relatively similar physi- 
cal parameters for most cores. We compare their results to the 
parameters obtained for the cores from our sample and test to 
which extent these results reflect the constraints provided by 
the set of lines observed, the limitations of the data analysis or 
the real physical properties of the cores. To study the influence 
of the data analysis we compare the escape probabiUty approxi- 
mation with a self-consistent radiative transfer code computing 
the excitation conditions in an inhomogeneous core with inter- 
nal turbulence. 

In Sect. 2 we provide a short overview on the sample and 
the observations. Sect. 3 discusses the traditional way to derive 
the cloud parameters from the observations. Using the mathe- 
matical description of the radiative transfer problem in Appx. 
A we compute cloud parameters in an escape probability ap- 
proximation. In Sect. 4 the fully self-consistent radiative trans- 
fer code from Appx. B is used to derive the cloud parameters. 
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Table 1. Selected source sample with central positions 



2.2. Observational details 



Source 




Q1950 




<5l950 




W49A 


19'' 


07"' 


50.0^ 


09° 


01' 


32.0" 


W33 


18'' 


11™ 


19.8^ 


-17° 


56' 


21.0" 


W51A 


19'' 


21"' 


25.4^ 


14° 


24' 


45.0" 


W3(0H) 


2h 


23™ 


15.4^ 


61° 


38' 


53.0" 


W3 


2h 


21™ 


42.5^ 


61° 


52' 


22.0" 


S255 


6" 


09™ 


58.4^ 


18° 


00' 


19.0" 


S235B 


5" 


37™ 


31.8^ 


35° 


40' 


18.0" 


S106 


20'' 


25™ 


39.8^ 


37° 


12' 


46.0" 


Serpens 


18'' 


27™ 


25.0^ 


01° 


12' 


00.0" 


DR21 


20'' 


37™ 


14.0^ 


42° 


09' 


00.0" 


MonR2 


6" 


05™ 


14.0^ 


-6° 


23' 


00.0" 


NGC2264 


6" 


38™ 


26.0^ 


9° 


32' 


00.0" 


OMC-2 


5" 


32™ 


57.0^ 


-5° 


12' 


12.0" 


pOph A 


16'' 


23™ 


17.9^ 


-24° 


17' 


18.0" 


NGC2024 


5" 


39™ 


12.0^ 


-1° 


56' 


40.0" 



Sect. 5 compares the resulting data with parameters obtained 
from independent observations and discusses implications for 
the physics of massive cores. 



2. Observations 

2.1. The sample 

The selection of the sample of massive cores was determined 
by the need of bright "standard" sources for the SWAS satel- 
lite which has a spatial resolution of about 4 arcmin in the fre- 
quency range between 490 and 560 GHz. For a comparison to 
data obtained at similar angular resolution we used observa- 
tions taken with the array receiver of the FCRAO 14m tele- 
scope at 98 GHz (about 1 arcmin resolution) and performed 
complementary observations with the 3m KOSMA telescope 
at 245 and 343 GHz where its beam size is approximately 2 ar- 
cmin. 

The sources were selected from the SWAS source list 
(Goldsmith et al., priv. comm.) to be observable from the 
FCRAO and KOSMA and bright enough to be detectable in a 
reasonable integration time. Tab. [l| lists the 15 sources selected 
with their central position. 

NGC2024 is not part of the SWAS sample but we have 
observed this region as a standard for comparison: it is rela- 
tively close (w 450 pc) and has been studied already by nu- 
merous authors using various techniques. An extended map 
of t he cloud and its environment in CS 2-1 was provided 
by Lada et al. (1991)t and line profiles in four transitions 
of the main CS iso tope are given by Lada et al. (1997). 
Mezger et al. (1992) have identified seven clumps in NGC2024 
from dust observations using the IRAM 30m telescope whereas 
the FCRAO and KOSMA beams can only distinguish between 
the two bright clumps FIR3 and FIR5. The data analysis is per- 
formed for the position of the brightest clump FIR5. Here, the 
KOSMA beam also contains weak contributions from FIR4, 
FIR6, and FIR7. 



All sources except NGC2024 were observed in CS 2-1 by 
Howe (priv. comm.) using the FCRAO 14m telescope providing 
a resolution of 53". and a main beam efficiency r/j^^^^ = 0.58. 
The cores were covered by 30-point maps with a sampling of 
50". The CS 2-1 spectra for the southern core in NGC2024 
were taken from Lada et al. (1997). 

The CS 5^ and 7-6 observations used the dual channel 
KOSMA SIS receiver with noise temperatures of about 95 K 
in the 230 GHz branch and 120 K in the 345 GHz branch. 
The 3m telescope provides a spatial resolution of 110" in CS 
5-4 and 80" in CS 7-6. The default observing mode for all 
cores were cross scans with a separation of 50" between sub- 
sequent points. Only for NGC2024 complete 5'x5' maps were 
obtained. At the time of the observations the telescope surface 
provided main beam efficiencies Ty^j^"* = 0.54 and rj"^^^^ = 
0.48 respectively for the two transitions. For all measurements 
we use the conservative estimate of about 10% uncertainty for 
the main beam efficiency, another 10% atmospheric calibration 
uncertainty and add another 5% for possible drifts etc. As sys- 
tematic errors they might sum up linearly to a total calibration 
error of at most 25%. 

The FCRAO spectrometer had a channel width of 19.5 kHz 
corresponding to a velocity spacing of 0.060 km/s. For the 
broad lines from the sample four velocity channels were 
binned. The resulting r.m.s. falls between 0.2 and 0.45 K. The 
KOSMA spectra were taken with the medium resolution spec- 
trometer (MRS) and the low resolution spectrometer (LRS) 
providing channel widths of 167 kHz and 688 kHz, respec- 
tively. Depending on the different combinations of these back- 
ends with the receivers at 245 and 343 GHz we obtain velocity 
spacings between 0.15 and 0.84 km/s. The particular spacing 
is not important for the analysis performed here because none 
of the lines shows strong spectral substructure. All points were 
integrated up to a noise limit of 0. 1 K per channel. 

2.3. Results 

For most sources the CS 2-1 maps show an approximately el- 
liptical intensity peak with a weak elongation at scales of a 
few times the resolution. In W49A, Serpens, DR21, Mon R2, 
OMC-2, p Oph A, and NGC2024 we can distinguish a second 
intensity maximum apart from the central position. Tab. ^ sum- 
marises the parameters of the line profiles at the central position 
for all cores. The majority of line profiles are approximately 
Gaussian as indicated by integrated line intensities close to the 
Gaussian value of 1.06 TmbAw in Tab. |[ Broad wings are only 
visible in S255, W33, and DR21. 

W49 shows a double-peak structure which has been inter- 
preted e.g. by pickel & Auer (1994) as the footprint of large 
scale collapse. They fitted HCO+ line profiles by a spherical 
collapse model but concluded that additional components are 
needed to explain the observations. Using the radiative transfer 
code from Appx. B we have tested their infall model and found 
that, while reproducing the HCO+ profiles, it completely fails 
to explain the CS observations. The enhanced blue emission 
characteristic for collapse is visible only in CS 7-6. In CS 2-1 
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Table 2. Measured line parameters at the central core positions 
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CS 2- 


-1 




CS 5^ 




CS 1~ 


-6 


Source 




J Ti^\hdv 


Av 


Ti-nh 




Av 




J Tyni^dv 


At; 




[K] 


[K km/s] 


[km/s] 


[K] 


[K km/s] 


[km/s] 


[K] 


[K km/s] 


[km/s] 




'ill 


35.17 


8.75 ± 0.55 


1.54 


12.61 


7.69 ±0.19 


1.23 


12.15 


9.28 ± 1.24 


W49A 


4.34 


32.76 


7 10 + 38 


2.04 


1 S 70 
1 o. / u 


8 67 + f) 1 ^ 


1.06 


9.52 


8 43 ± 1 03 


W33 


9.74 


67.24 


6 49 + 1 6 


3.76 






2.94 


17.44 


5 57 ± 05 


W51A 


10.86 


110.51 


Q 56 + 09 


4.99 


fi7 SQ 


1 7 7S+ n flQ 
iz,. / ozn U.U7 


3.38 


40.63 


11 78 ± OS 




5.50 


26.48 


4 5^ + 70 


2.13 


1 n 1 ^ 

iU. i J 


H-.HU ZIZ U.UO 


1.25 


5.81 


4 40 ± ^ 1 


W3 


1 .(fl 


38.97 


4 77 + 05 


1.35 


7.72 


5.36 ± 0.28 


1.56 


8.44 


5 08 ± 09 


S255 


8.02 


21.55 


2.53 ± 0.04 


2.91 


10.13 


3.28 ± 0.03 


1.67 


5.29 


2.99 ± 0.15 


S235B 


6.38 


16.21 


2.39 ± 0.05 


0.85 


2.67 


2.97 ± 0.13 


0.42 


1.10 


2.53 ±0.19 


S106 


3.29 


8.19 


2.34 ±0.10 


0.83 


2.11 


2.40 ± 0.22 


0.65 


2.23 


3.23 ± 0.78 


Serpens 


3.74 


9.36 


2.35 ±0.13 


0.67 


2.52 


3.59 ± 0.35 


0.39 


1.68 


4.11 ± 0.65 


DR21 


6.48 


25.00 


3.62 ±0.10 


2.41 


9.74 


3.80 ±0.18 


2.92 


10.92 


3.51 ±0.06 


Mon R2 


5.53 


12.26 


1.90 ±0.16 


2.15 


6.52 


2.38 ± 0.27 


1.70 


3.83 


1.95 ±0.19 


NGC2264 


5.64 


21.90 


3.66 ± 0.08 


2.30 


9.48 


3.89 ± 0.08 


1.25 


5.42 


3.91 ± 0.12 


OMC-2 


5.62 


9.26 


1.55 ± 0.05 


1.96 


5.93 


2.85 ±0.10 


0.92 


1.54 


1.58 ±0.25 


pOph A 


4.22 


11.10 


2.47 ± 0.25 


2.07 


2.30 


1.04 ±0.08 


<0.2 







Table 3. Measured line parameters within NGC2024 



Transition 


[K] 


J Tmhdv 
[K km/s] 


Av 
[km/s] 


C^S 5^ 


5.19 


11.65 


2.10 ±0.02 


C^^S 7-6 


3.63 


9.13 


2.38 ± 0.08 


C^^S 5-4 


1.19 


2.22 


1.75 ±0.09 


C^^S 7-6 


0.65 


0.98 


1.45 ± 0.13 



and 5^ we rather find an enhanced red emission. No spheri- 
cally symmetric collapse model can explain these observations. 
Instead of constructing a more complex model we have simply 
decomposed the emission into two separate components with 
a relative velocity of 8.5 km/s in the line of sight denoted as 
W49A('') and W49A(''). From the modelling in Sect. 4 it turns 
out that we cannot even distinguish whether the two compo- 
nents are moving towards each other or apart as long as they do 
not line up exactly along the line of sight. Hence, we will treat 
them separately in the following, ignoring any possible inter- 
action. Further observations including other tracers should be 
included to better resolve the situation. 

In NGC2024 we also mapped the less abundant isotope 
C^'^S with KOSMA in addition to the main CS isotope ob- 
served in all cores. The observed line parameters are given in 
Tab. ^ Lada et al. (1997) provided detailed CS spectra for the 
southern core at the position of the FIRS. They obtained in the 
CS 2-1 transition T,„b = 15.4 K, Au = 1.80 km/s at 24" 
resolution, in CS 5^ T,nb = 9.6 K, Aw = 2.1 km/s at 30" 
resolution, in CS 7-6 T,nb = 10.1 K, Av = 2.7 km/s at 20" 
resolution, and in CS 10-9 T,nb = 10.6 K, Av = 2.1 km/s at 
14" resolution. 



3. Cloud parameters from the escape probability 
model 

To interpret the small amount of information contained in ob- 
servations of at most five transitions showing mainly Gaussian 



profiles and essentially unresolved approximately circular sym- 
metric intensity distributions, we need a simple cloud model 
that is both physically reasonable and characterised by few pa- 
rameters. An obvious choice is a spherically symmetric model. 
This geometry reflects early phases and t he large scale be - 
haviourof several collapse simulations (e.g. Galli et al. 1999 ), 
whereas the inner pa rts of collapsing clouds are probably flat- 
tened structures (e.g. Li & Shu 1997 ). 

Even in spherical geometry there is no simple way to solve 
the radiative transfer problem relating the cloud parameters to 
the emitted line intensities (see Appx. A.l). Thus we cannot 
compute the cloud properties directly from the observations. 

3. 1. Application ofttie escape probability 
approximation 

A common approach is the escape probability approximation 
discussed in detail in Appx. A. 2. Assuming that all cloud pa- 
rameters, including the excitation temperatures, are constant 
within a spherical cloud volume one can derive a simple for- 
malism relating the three parameters kinetic temperature Tkin, 
gas density n-a^ , and column density of radiating molecules on 
the scale of the global velocity variation A^moi/Ati to the line 
intensity at the cloud model surface. No assumption on molec- 
ular abundances is required. 

Since a telescope does not provide a simple pencil beam 
we have to correct the model surface brightness temperature 
by the beam filling factor rjf, given as the convolution inte- 
gral of the normalised intensity distribution with the telescope 
beam pattern, to compute the observable beam temperature. 
Unfortunately, the brightness profile of the source is a non- 
analytic function where we can only give simple expressions 
for the central value observed in a beam much smaller than the 
source or for the integr al va lue observ ed in a beam much larger 
than the source (Eqs. and \A.ll\ ). For intermediate situa- 
tions we approximate the beam temperature by starting from 
both limits and using a beam filling factor given by the con- 
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Table 4. Source size corrected for beam convolution. 



Source 


FWHM in a 


['] 


FWHM in S 


['] 




2-1 


5-4 


7-6 


2- 


1 5-4 


7-6 


W49A 


1.0 


1.1 


0.3 


1.5 


1.7 


1.3 


W33 


1.7 


1.7 


1.4 


1.7 


1.3 


1.0 


W51A 


1.9 


1.7 


1.7 


2.1 


2.2 


1.1 


W3(OH) 


2.1 


1.1 


1.3 


1.7 


1.4 


1.4 


W3 


1.0 


0.6 


0.8 


1.0 


0.6 


0.3 


S255 


1.9 


1.6 


1.1 


1.2 


0.6 


1.1 


S235B 


1.0 


2.2 


1.1 


0.9 


0.2 


0.3 


S106 


3.0 


1.6 




2.5 


2.2 




Serpens 


2.0 


2.2 




2.1 


2.3 




DR21 


3.0 


2.6 


1.8 


1.5 


1.7 


1.8 


Mon R2 


3.2 


3.0 


1.1 


3.6 


3.2 


1.6 


NGC2264 


3.2 


1.7 


0.6 


2.9 


2.1 




OMC-2 


2.4 


1.4 


0.6 


2.5 


1.4 


0.6 


pOphA 


1.3 


3.0 


2.6 


1.9 


2.9 


2.5 


NGC2024 


1.1 


1.1 


1.3 


1.3 


2.5 


2.3 



volution integral of two Gaussians. The difference between the 
two values provides an estimate of the error made in the beam 
convolution. 

To compute the integral we fitted the observed brightness 
distributions by Gaussians. Most cores are well approximated 
by slightly elongated Gaussians. W49A, S235B, Serpens, and 
p Oph A show asymmetric scans so that the size determina- 
tion is somewhat uncertain. The fit error is about 0.3' for these 
three sources. For the rest of the cores we obtain typical val- 
ues of less than 0.2'. The true object size finally follows from 
the deconvolution of the measured intensity distribution with 
the telescope beam. The resulting source sizes in a and 6 are 
given in Tab. ^. As the geometric mean is sufficient to com- 
pute the beam filling factor we don't expect any serious error 
from the fact that the cross-scans in a and S do not necessarily 
trace the major axes of the brightness distribution. For sources 
which are considerably smaller than the beam widths of 53", 
107", and 80", respectively, only a rough size estimate is pos- 
sible according to the nonlinearity of the deconvolution. This 
holds for W49A, W3, S235B, and partially S255. Most clouds, 
however, show an extent of the emission which is close to the 
beam size. 

In general different values are obtained for the spatial 
FWHMs in the different lines. In Tab. Q we find two classes 
of sources with respect to the variation of the source size de- 
pending on the transition observed. Most cores show a mono- 
tonic decrease of the visible size when going to higher transi- 
tions. This is expected from the picture that higher transitions 
are only excited in denser and smaller regions. Serpens, p Oph 
A, and NGC2024, however show the smallest width of the fit 
in the CS 2-1 transition. This is explained by eye inspecting 
the 2-1 maps and corresponding high-resolution observations 
from the literature where we see that the three sources break up 
into several clumps which are only separated in the 53" beam 
but unresolved in the KOSMA beams. In these cases, we have 
restricted the analysis to the major core seen in the CS 2-1 
maps using its size to compute the beam filling, although this 
approach introduces a small error in the data analysis by as- 



signing the whole flux measured in the higher transitions to 
this central core. 

Applying the two limits for the beam size treatment (Eqs. 
A. 9 and |A.11[ ) using the sizes from Tab. ^ we find that the 
resulting gas and column densities are the same within 20% 
except for W49A, W3, S235B, and NGC2264. The first three 
are small compared to the beams so that the results from Eq. 
( |A.9| ) have to be rejected and only Eq. ( |A.11[ ) can be used. For 
NGC2264 we cannot provide a simple explanation for the dif- 
ference so that we give a relatively large error bar covering the 
results from both approximations. 

The size of the resulting parameter range in Tkm, nn^, and 
Nqs/^v is determined by the accuracy of the observations. 
For two cores it was only possible to set a lower limit to the gas 
density. Moreover, we were not able to provide any good con- 
straint to the cloud temperature for all sources. Values between 
about 30 K and 150 K are possible. Hence, an independent de- 
termination of the cloud temperatures is required. Several dif- 
ferent methods based on optically thick CO, NH3 or dust ob- 
servations are discussed in the literature and we used the values 
from the references given in Tab. |^. In addition to the se values 
we also used 50 K as assumed by Plume et al. (1997) as "stan- 
dard" temperature in the parameter determination for massive 



3.2. Resulting core parameters 

Tab. H lists the parameters from the escape probability model 
for all cores. Whereas the column density is well constrained 
for most clouds, there is a considerable uncertainty in the gas 
density resulting from the unknown cloud temperature. At the 
temperature of 50 K we obtain average values and logarithmic 
standard deviation factors of 



(rinj = 7.9 10^ 
(iVcs/Aw) = 1.21014 



x/^ 1.5 
x/~2.8 



(1) 



From the 71 cores analysed by Plume et al. (1997) 
ing this temperature they obtained 



(iVcs/Aw) 



8.5 10^ 
2.5 IQi-* 



x/~ 1.7 
x/-f 3.1 



(2) 



This good agreement indicates first that both investigations 
study the same type of clouds visible in the CS transitions. 
Second, this shows that the static escape probability (Eq. |A.l l[ ) 
used here and the LVG escape probability (Eq. |A.9[ ) applied 
by Plume et al. (1997) differ only marginally as discussed al- 
ready by Stutzki & Winnewisser (1985) Third, the observa- 
tional data from both telescopes give almost equivalent results, 
i.e. the reliability of the parameters hardly profits from using 
the 10 times smaller beam of the IRAM telescope when the 
escape probability approximation is used. 

For NGC2024 we are able to test the consistency of the 
results from the CS and the C'^^S observations. The result- 
ing hydrogen densities agree for both isotopes within 20% 
at all temperatures assumed. We obtain 9 10^ cm^^ at 25 K 
(Mezger et al. 1992) and 3.2 10*^ cm^^ at 40 K (Ho et al. 1993). 
Unfortunately, this is a core where we can only give a lower 
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Table 5. Clump parameters derived from the escape probability 
model 



Source 




"H2 




Ncs/Av 
[cm^^/kms-'] 


x/- 




[K] 


[cm-''] 








20" 


> 7.3 10** 




2.2 10''' 


1.3 




50" 


1.310'' 


1.5 


2.1 10^'' 


1.2 




20" 


4.5 10'^ 


1.6 


2.8 10''' 


1.3 




50" 


7.6 lO'^ 


1.5 


2.7 10'" 


1.2 


W33 


40'" 


1.7 lO'' 


1.7 


2.3 lO'"* 


1.2 




50"^ 


1.2 lO'' 


1.6 


2.3 10'" 


1.2 


W51A 


20" 


2.1 10^ 


1.6 


4.8 10'" 


1.7 




50 


1.1 lO'' 


1.5 


4.1 10'" 


1.4 




57" 


8.9 lO'^ 


1.4 


3.7 10'" 


1.4 


W3(OH) 


30'' 


1.9 lO'' 


1.5 


1.3 10'" 


1.2 




50 


8.5 10^ 


1.4 


1.3 10'" 


1.2 


W3 


30== 


7.1 lO'' 


2.9 


7.4 10'" 


1.5 




50 


9.8 lO'^ 


1.8 


8.1 10'" 


1.4 




55'' 


7.4 10^ 


1.6 


8.3 10'" 


1.3 


S255 


40*5 


1.3 lO'' 


1.6 


1.310'" 


1.3 




50 


9.3 10^ 


1.5 


1.210'" 


1.3 


S235B 


40'^ 


1.2 10'^ 


1.5 


3.7 10'^ 


1.3 




50 


8.9 lO'^ 


1.5 


3.6 10'=* 


1.3 


S106 


10' 


>6.3 10'' 




4.7 10'^ 


1.4 




25' 


7.1 10^ 


1.5 


3.8 10'^ 


1.4 




50 


2.5 10^ 


1.6 


3.910'^ 


1.4 


Serpens 


25^ 


4.6 10^ 


1.2 


2.8 10'^ 


1.5 




50 


1.910^ 


1.6 


2.8 10'^ 


1.5 


DR21 


35"^ 


1.5 10® 


1.7 


9.3 10'^ 


1.9 




50 


9.1 10^ 


1.6 


8.7 10'^ 


1.4 


Mon R2 


25' 


1.7 lO" 


1.8 


4.5 10'^ 


1.5 




50"" 


6.2 10^ 


1.7 


3.7 10'^ 


1.7 
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1.6 
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50 
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1.7 
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1.4 
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1.6 
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50 
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1.8 
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NGC2024 


25" 
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40' 
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50 
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1.7 
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Sievers et al. (1991) , ^ Goldsmith & Mao (1983) 



Haschick & Ho (1983), '' Wilson et al. (1991) 



rieftrunk et al. (1998]|, ^ [Tieftrunk et al. (1995) 
Faffe et al. (1984j, [Nakano & Yoshida (1986) 



Roberts et al. (1997; , ^ McMulhn et al. (1994) 



Garden & Carlstrom 1992, ' 


Montalban et al. (1990;, 


™ iGiannakopoulou et al. (1997), " 


Krugel etal. (1987) 



1 



Castets & Langer (1995)|, ^ patrla et al. (1983)|. 
'I ICesaroni & Wilson (1994)1 ' Liseau et al. (1995 



Ho et al. (1993; , ' Mezger et al. (1992) 



limit to the column densities. The limits deviate by a factor 13, 
which is significantly different from the ter restrial isotopic ratio 
of 23 but close to the value 10 derived by Mundy et al. (1986) 
for the isotopic ratio in NGC2024. 

The escape probability model provides a first estimate to 
the physical parameters but its limitations are obvious. It is 



definitely not justified to assume constant parameters within 
the whole cloud. Moreover, several observations are in contra- 
diction to the parameters from the escape probability models. 
pada et al. (1997)| detect ed the CS 10-9 transition in NGC2024 
and Plume et al. (1997) observed the 10-9 and 14-13 transi- 



tions in S255 and W3(OH). The critical densities for these 
transitions are about 6 10^ cm^^^ and 2 lO'* cm^^^ respectively. 
From the densities in Tab. ^ one would conclude that these 
transitions are not excited. Hence, a more sophisticated model 
has to be applied to obt ain a physically reas onable explana- 
tion of the measurements. Plume et al. (1997) suggested a two- 
component model or continuous density gradients to resolve 
this contradiction. We will discuss a self-consistent radiative 
transfer model including a radial density profile in the follow- 
ing. 



4. Cloud parameters from the nonlocal model 

4. 1. Line fitting by SimLine 

We performed nonlocal radiative transfer simulations using 
the line radiative transfer code SimLine introduced in detail in 
Appx. B. SimLine is a FORTRAN code to compute the profiles 
of molecular rotational lines in spherically symmetric clouds 
with arbitrary density, temperature and velocity distribution. It 
consists of two parts: the self-consistent solution of the balance 
equations for all level populations and energy densities at all 
radial points and the computation of the emergent line profiles 
observed by a telescope with finite beam width and arbitrary 
offset. The optical depths in the lines may vary from minus a 
few, corresponding to weak masing, to several thousand. 

Already in the spherically symmetric description of a core 
we face a large number of parameters. For all quantities (hydro- 
gen density, kinetic temperature, velocity dispersion, molecular 
abundances) a radial function has to be found. Regarding the 
limited amount of information available from the three transi- 
tions, this leaves many options open. We decided to assume 
simple power-law radial functions and a central region with 
constant parameters for all quantities in the core simulations. 
This reduces the number of parameters to two (central value 
and radial exponent) for each gas property, plus the outer and 
inner radius. 

The parameter fit procedure used the multidimensional 
downhill simplex algorithm from Press et al. (1992). Although 
it is not the most efficient way in terms of convergence speed 
it turned out to be very robust in all situations considered. 
Because a downhill simplex code does not necessarily find the 
global minimum of the function we performed for each core 
several runs with randomly chosen initial simplex covering a 
large part of the physically reasonable parameter space. For all 
clouds we made at least 30 runs to get a rough idea of the topol- 
ogy of the function. For cores like W33 this turned out to be 
sufficient since only one large minimum showed up which was 
found in almost half of the runs. The other extreme is S106 
where we needed almost 1000 runs to be sure that we found 
the global minimum. Here, the function was quite complex 
with numerous local minima. Future improvements of the fit 
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procedure should include more sophisticated algorithms like 
simulated-annealing approaches. 

The noise in the line profiles produces some graininess of 
the function when directly fitting the measured profiles. This 
results in a very slow convergence of the algorithm close to the 
minimum. A considerable acceleration can be obtained by not 
fitting the measured noisy data but a smooth approximation to 
them. The line profiles were represented by a superposition of 
a Gaussian and a Lorentzian profile. This allows a good char- 
acterisation of all measured profiles including the reproduction 
of asymmetric profiles, self-absorption dips and line wings. 

Taking the three measured transitions and their spatial ex- 
tent as the quantities to be reproduced by a fit we find that 
we are able to derive at most six parameters to a reasonable ac- 
curacy. Fits with seven or eight free parameters, although still 
slightly improving the numerical value, do not produce any 
significant changes above the noise limit. This is comparable 
to the results by Young et al. (1998) fitting the full position ve- 
locity map of a particular core in a single transition. There is 
however a number of additional parameters where the model 
can derive certain limits. Hence, we have to decide first which 
parameters should be fitted directly, which parameters may be 
constrained by preventing reasonable fits when outside a cer- 
tain range, and which parameters can be guessed independently 
from a physical line of reasoning. This is discussed in detail for 
all quantities in the following subsections. 

Unfortunately, it is impossible to give an comprehensive 
error estimate for the parameters derived from the simulations. 
This would need a description of the six-dimensional surface 
in the parameter space within which none of the observational 
error bars is exceeded. According to the complex topology of 
the function it is not possible to give an easy description for 
the six-dimensional valley around the global minimum or its 
boundaries. 

As a simple alternative we performed only one-dimensional 
variations to get a rough estimate of the maximum error in the 
parameters that we must expect. After the convergence of the 

fit, each of the fit parameters was varied up and down un- 
til one of the computed lines deviated by the assumed maxi- 
mum observational error of 25 % from the measured lines. The 
central values of the different functions were varied indepen- 
dently. When changing the inner radius, the central values of 
the gas parameters were adjusted to keep the functions in the 
power-law region unchanged. When varying exponents the cor- 
responding central values were corrected in such a way that the 
parameters at the density of 2 lO^cm^'^ remained constant. By 
the selection of this density as the fix point when changing the 
slope, we scan about the maximum possible range of the expo- 
nents. This provides a conservative estimate of the maximum 
error 



4.2. The turbulence description 

In spite of the large number of free parameters in the cloud 
models it is impossible to fit the line profiles with a smooth den- 



case. Thus we have to take into account the effects of internal 
clumping and turbulence leading to a more realistic picture and 
preventing strong self-absorption. 

SimLine treats turbulence and clumping in a local statis- 
tical approximation following Martin et al. (1984) (see Appx. 
B.2). The cloud material is subdivided into small coherent units 
(clumps) with only thermal internal velocity dispersion. The 
relative motion of many units then provides the full velocity 
profile. Martin et al. (1984) showed that the effective optical 
depth of such an ensemble can be considerably reduced com- 
pared to the microturbulent approximation. The relative reduc- 
tion of the total optical depth depends on the optical depth, i.e. 
of size the coherent units. 

This description does not need any assumptions on the na- 
ture of the turbulence creating the internal cloud structure. The 
reduction occurs in the same way whether clumps are units 
of the same velocity in a medium of constant density, repre- 
senting the behaviour of incompressible turbulence, or whether 
they are density enhancements in a thin inter-clump medium. 
The different nature of these two scenarios has to be taken into 
account, however, when computing the excitation. Both the av- 
erage gas density, providing a measure for the column density 
and thus the line intensity, and the local density, providing the 
coUisional coupling to the gas, enter the balance equations. In 
case of coherent units in velocity space both densities agree. 

For density clumps we use the additional simplification to 
treat the cloud locally as a two-component medium neglect- 
ing the contribution of the inter-clump medium to the radiative 
transfer Then the collisional excitation is provided by the den- 
sity in the clumps and the column density is provided by the av- 
erage density. The ratio between the two quantities reflects the 
filling factor of the volume occupied by dense clumps. From 
the mathematical point of view this is equivalent to the treat- 
ment of the cloud as a homogeneous medium where the clump- 
ing occurs only in velocity space and the abundance of the ra- 
diating molecules is reduced. Hence, it is impossible to sepa- 
rate the influence of the filling factor from that of the molecular 
abundance so that the line fit provides only a combined quantity 
which we will denote combined abundance in the following. 
Beside the modification of the molecular abundance the statis- 
tical turbulence description introduces the size of the coherent 
units as an additional parameter 

Fig. |l] demonstrates the influence of the turbulence. It 
shows the central line profiles of the best fits to the S255 
observations using different assumptions on the turbulent na- 
ture. We selected S255 here because we can exploit the ad- 
va ntage of additional d ata for the CS 10-9 transition measured 
by Plume et al. (1997) providing more constraints on the cloud 
model[[ In all models the core is optically thin in CS 7-6 and 
10-9, so that the reduction of the optical depth by the turbulent 
clumping produces only minor changes in these lines. In CS 
2-1 the differences are, however, most obvious. In the micro- 
turbulent description and for the incompressible turbulence we 
find self-absorbed line profiles. A reasonable fit to the obser- 



sity and velocity distribution. The ty pical self-ab sorbed profiles 
known for all microturbulent codes ( [White 1977 ) appear in this 



' Unfortunately, [Plume et al. (1997) listed only the peak intensity 
and FWHM so that we have to assume a Gaussian for the CS 10-9 
profile here. 
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Fig. 1. Line profiles observed at the central position of S255 in 
four CS transitions together with best-fit models. The dashed 
line represents the microturbulent approximation, the dotted 
line stands for clumping in velocity space with cell sizes of 
0.01 pc and the solid line includes the additional effect of den- 
sity clumping where the combined abundance is reduced by a 
factor 100. 



vations is only possible using the turbulence model including 
clumps with enhanced density. The reduction of the effective 
optical depth of the cloud by increasing the optical depth of the 
coherent units produces narrower lines providing a better fit to 
the observed line profiles. 

We find that the parameter fits do not provide accurate val- 
ues for the clump size and the combined abundances reflecting 
the volume filling factor. They only constrain an interval of pos- 
sible values. We obtained good fits to the observations for the 
full range of clump sizes between about 0.005 and 0.05 pc. This 
size scale is confirmed by independent determinations of clump 
sizes by high-resolution observations in some of our cores. For 



NGC2024 Mezgeret al. (1992) determined a radius of dense 



conde nsations around 0.015 pc and the interferometric stud- 
ies by Wiesemeyer et al. (1997) showed val ues between 0.005 
and 0.01 pc. In W3 Fieftrunket al. (1998) observed compact 



clumps with a size of 0.02 pc and in OMC-2 |Chini et al. (1997)| 



found dust condensations with radii between 0.01 and 0.05 pc. 

For a better comparison to turbulence theory we prefer to 
specify the cell size in terms of the turbulent coiTelation length 
Icon which should be on the order of 0. 1 pc (Miesch et al. 1994, 
Goodman et al. 1998). The size of the units which are coher- 
ent with respect to the line radiative transfer is smaller by the 
ratio of the thermal line width to the total velocity dispersion. 
Thus the size range found corresponds to correlation lengths 
between 0.04 and 0.4 pc. In the following computations we use 
the intermediate value of 0.1 pc as correlation length for all 
clouds. 

Regarding the combined abundances we find two classes 
of objects. The majority of cores, including the example of 
S255, is weU fitted by values between 10~" and 10"^", 
whereas a second class, consisting of W3, Serpens, and S106 
needs values between 10"^'^ and 10~^ for a good fit. Because 
we do not know the molecular abundances there is no way 



to translate these values directly into clump filling factors. 
Assuming the CS abundances of 1.3 10"^ to 1.3 10^® obtained 
by Hatchell et al. (1998) for several star-forming cores, the first 
class corresponds to filling factors around 0.01, whereas the CS 
abundance of 4 10^^° from Plume et al. (1997) corresponds to 
a filling factor of 0. 1 . The filling factors in the second group are 
ten times higher accordingly. In all following computations we 
have used a combined abundance factor of 3 10^^^ for the first 
and 3 10^^° for the second group. In the translation to cloud 
masses we will use the intermediate CS abundance of 4 10~^. 

As an additional parameter quantifying turbulence in a 
one-dimensional cloud model we have to take a radial vari- 
ation of the turbulent velocity distribution into account to 
explai n the observ ed size-line width and size-density rela- 
tions ( Larson 1969 ). Exponents of this radial dependence be- 
tween about 0.1 and 0.7 are typically discussed (see e.g. 
[Goodman et al. 19^ ). We left the width and the exponent of 
the turbulent velocity distribution as free parameters and ob- 
tained exponents between 0.15 and 0.65. 



4.3. The density structure 

Any physically reasonable cloud model should include a 
spatial dependence of the gas parameters. Collapse sim- 
ulations might provide reasonable model assumptions for 
the density structure. Bodenheimer & Sweigart (1968) and 
[Shu (1977) have shown that an isothermal sphere evolves into 
a power law density profile nn^ oc r^"" with Q!„ — 2.0. 
Homologous collapse simulations provide an exponent q;„ = 
3.0 (Dickel & Auer 1994) and the free-fall collapse discussed 



by Welch et al. (1987) results in a density exponent q:„ = 1.5. 



The inside-out collapse model (Shu 1997) combined two re- 
gions of different exponents and recent more sophisticated col- 



lapse simulations (see e.g. Basu & Mouschovias 1995) show 



more complex density structures with an average exponent a„ 
between 1.5 and 1.7. Dust observations of the density profile 
of protostellar cores show evidence both for cores with a typ- 
ical profile and for cores with fla t density structure and 



a sharp outer edge ( Andre et al. 1999 ). Thus we expect ex- 
ponents between about 1.5 and 2.0 in our power-law density 
model which is bound by an outer cut-off and a central con- 
stant region. Although this simple model may not reflect the 
whole complexity of the density profile we can hardly derive 
any more information from the limited observations available. 
We left the central density and the density exponent as free 
parameters to be fitted. The resulting exponents span the rela- 
tively wide range between 1.1 and 2.2. 

In the fit of the radii confining the power-law density pro- 
file we face two problems. We cannot distinguish changes of 
the model parameters on the smallest size scales where even 
the highest observed transition is thermalised because of the 
high density. Thus we can only set an upper limit R"l^^ to the 
radius of the central region where a transition from the power 
law behaviour in the envelope to constant parameters might oc- 
cur. Only in four clouds - Serpens, Mon R2, p Oph A, and the 
b-component of W49A - the strength of the CS 7-6 fine sets an 
upper limit to the density so that we can derive the inner radius 
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Table 6. Resulting cloud parameters for S255 assuming differ- 
ent exponents for the radial temperature dependence and fitting 
only the central Une profiles. 







(T) 






-^cloud 


Afcloud 


FWHM^ 




[K] 


[K] 




[cm--^] 


[pc] 


[Mq] 


['] 


0.0 


47 


47 


-1.5 


3300 


1.1 
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1.8 


-0.12 


73 


52 


-1.3 


2500 


1.3 


1500 


2.2 


-0.4 


220 


52 


-1.0 


1000 


2.5 


7000 


3.6 



CS 2-1 transition 

directly from the observations. In all parameter fits, we left the 
central radius as a free parameter. After the fit we increased the 
radius until the maximum deviation in one of the lines reached 
5 %. This provides a reasonable upper limit to the inner radius 
in all cases where the line profiles are independent of the den- 
sity structure below this limit and gives only small modifica- 
tions in the four cases where the inner radius was already well 
constrained by the fit. 

The outer radius of the cloud is also quite uncertain. We can 
easily provide a value for the extent of gas at densities above 
about 3 10^ cm^"^ based on the spatial extent of the CS 2-1 
emission and the line profiles. But it is not possible to derive 
a reliable value for the extent of low density gas. We can only 
give a lower limit to the outer radius and thus the mass of the 
massive cores considered. As the density exponents are typi- 
cally shallower than -2, a majority of mass could be present 
beyond this radial limit at low densities invisible in CS. We 
excluded the outer radius from the parameter fit using a suffi- 
ciently large value for all clouds and performed later a separate 
run reducing the radius to find the minimum outer radius i?™" ^ 
in a way equivalent to the maximisation of the inner radius de- 
scribed above. 

4.4. The temperature structure 

The temperature distribution of massive cores is still a mat- 
ter of debate (cf. Garay & Lizano 1999| ). During early phases 
of cloud collapse the temperature should remain constant as 
long as the core remains optically thin. Deviations are to be ex- 
pected, however, as protostellar sources are formed in most of 
our cores, leading to an internal heating of the cloud. Moreover 
in thin outer regions external heating can be important. Based 
on several observational results Scoville & Kwan (1976) set up 
a spherical cloud model with a warm inner region resulting in 
an temperature profile Tkin oc r"^ with ar ~ —0.4. 

Hence, we should also derive the core temperature and the 
temperature exponent from the radiative transfer model. In a 
first run we have investigated the influence of a temperature 
gradient in S255 when fitting only the central line profiles. We 
compared the best fit models to the S255 observations using ei- 
ther a constant cloud temperature, a temperature decaying with 
the exponent -0.4, or the temperature exponent as a free param- 
eter. The latter case provided a best fitting exponent of -0.12. 
All three fits showed an excellent agreement in terms of the 
central line profiles falling almost exactly on the solid curves 
in Fig. |l]. Thus, it is impossible to favour a certain exponent 
from the least squares fit of the line profiles only. 



In Tab. |g we see the resulting values for the other core 
parameters. The average temperature given here is computed 
as the mass-weighted average up to i?™oud- main differ- 
ence we find a kind of compensation between temperature 
and density exponent. The sum of both exponents is kept 
approximately constant to fulfil the constraints given by the 
line ratios. Temperature and density gradient act in a similar 
manner, leading to higher excitation in regions which are ei- 
ther denser or warmer. The variation of the density gradient, 
however, changes the extent of the emission in the CS 2-1 
transition (the change is much smaller in the higher transi- 
tions). Consequently, we can constrain the temperature expo- 
nent when taking the observed size of the source into account. 
The beam convolved FWHM of the CS 2-1 emission in S255 
falls between 1.5' and 2.1' clearly excluding models with steep 
temperature gradients. 

In the fit procedure applied to all cores we have thus in- 
cluded the fit of the spatial extent. In this way we can constrain 
the temperature exponent but we are not able to derive exact 
values. In the parameter fit we always start with an isothermal 
cloud and change the temperature gradient in steps of 0. 1 un- 
til the model provides a simultaneous good fit to the central 
line profiles and the spatial FWHM of the emission. It turns 
out that the observations of all clouds except W33, W3 and 
Mon R2 can be fitted by an isothermal model. This indicates 
that a large part of the gas mass in the clouds is characterised 
by a uniform kinetic temperature. This isothermal behaviour, 
however, does not support the assumption of constant excita- 
tion temperatures like in escape probability model. The excita- 
tion temperatures show local variations which are steeper than 
the density gradient. In the derivation of the kinetic tempera- 
ture structure we have to keep in mind, however, that we are 
not very sensitive to the exact value of the kinetic tempera- 
ture gradient. Observations of higher transitions like CS 14-13 
would be needed to obtain reliable values including full error 
estimates. 



4.5. Other parameters 

Beside the temperature gradient we can expect the formation 
of compact HII regions in the centre of a star forming core. 
Using the parameters of the central HII region derived by 
|Dickel&Auer( 1994)1 for W49 (i? = 0.2 pc, = 2.6 10^ 



lO** K) we have compared the resulting CS lines 
when either including or neglecting the HII region in the ra- 
diative transfer computations (see Appx. B.3). We find that the 
influence of the HII region is negligible in this example for the 
four CS transitions considered here. 

In general HII regions have two dominant effects on the 
line profiles. The integrated brightness changes the molecu- 
lar excitation throughout the cloud and is visible as contin- 
uum emission. Moreover, the molecular material in front of an 
HII region appears partially in absorption. A large HII region 
with electron densities being a factor 10 or more higher than 
in the example above, would thus result in distinct changes in 
the lines. Depending on the configuration and velocity struc- 
ture, the CS lines may appear in absorption or with P Cygni 
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Fig. 2. Observed line profiles and best fit model to the obser- 
vations of NGC2024 using a CS to C^^'S abundance ratio of 
13. 

profiles. Moreover, a strong continuum emission would be ob- 
served at all frequencies considered. However, bright HII re- 
gions w ith high electron densities are unlikely to be that ex- 
tended ( Wood & Churchwell 1989 ). For ultracompact HII re- 
gions the change of the line profiles by absorption is negli- 
gible due to the small angular size of the HII region so that 
the molecular excitation is the remaining effect. From the lack 
of a bright continuum underlying the lines we can, however, 
exclude configurations with a bright HII region here. Weaker 
compact or ultracompact HII regions - although possible - 
were not included in the fit computations as they would only 
influence material in their close environment which cannot be 
resolved in this study. 

Regarding the chemical evolution of massive cores one 
should also expect a variation of the molecular abundances of 
CS and C^^S. However, our present knowledge is still insuffi- 
cient to guess reliable values here (cf. Bergin & Langer 1997). 
As discussed in Sect. [l|, it is also not possible to fit the abun- 
dance independently from the clump filling factor, so that we 
adopted here a constant CS abundance of 4 10^^ relative to H2 
ignoring any radial variation of the abundances. 

An additional test is only possible for NGC2024 where 
C^^S is sufficiently bright so that we could include it in the 
fits. Fig. ^ shows the best fit model to all six available line 
profiles. Equivalent to the results from the escape probability 
model we get the best match for a relative molecular abundance 
X(CS)/X(C^''S) w 13. AU six lines are simultaneously fitted. 

None of the cores except W49A show clearly asymmet- 
ric profiles as a signature of systematic internal velocities. The 



weak wings seen in a few of the other lines at low resolution 
are insufficient to derive any collapse or outflow model. Thus 
we have fitted all cores with a static model. 



4.6. Comparing the massive cores 

Tab. ^ lists the resulting best fit parameters for all cores with 
their eiTor intervals. As discussed in Sect. 4.3 the intervals for 



the radii are in most cases only limited at one end. The central 
density refers to the value at the radius i?™™^- 

Comparing the different cores we find that the inner and 
outer radii are mainly determined by the selectional bias from 
the observability as one massive core. The most distant cores 
are only detectable with the KOSMA telescope if they are rel- 
atively large whereas at small distances only small cores are 
unresolved. 

The minimum central density fitting the CS lines varies be- 
tween 10^ cm-3 for W49A('') and 8 10^ cm^^ for NGC2024. 
The high value derived for NGC2024 is due to the availabil- 
ity of the CS 10-9 observations tracing higher densities. With 
the three CS lines measured for most cores, only the density 
range below about 10^ can be reliably traced, so that we have to 
take the central density for all clouds except W49A(''\ Serpens, 
Mon R2, and p Oph A as lower limits. The density exponent 
covers the range between -1.1 and -2.2 where the majority of 
clouds shows values around -1.6 corresponding to large-scale 



collapse models (Sect. 4.3). 

Although we don't have a sample where we can expect 
to set up statistically significant correlations we can interpret 
some general relations. It turns out that the parameters are not 
completely independent of each other. Whereas the majority of 
clouds shows an average temperature between 20 and 50 K, the 
three clouds which needed a larger combined abundance fac- 
tor corresponding to a higher clump filling factor in the turbu- 
lence description (W3, S106, Serpens) also tend to require rel- 
atively high temperatures. There is a clear correlation between 
the cloud temperature and the turbulent line width indicating 
that heating and turbulent driving might have a related cause. 
The cores with a significant temperature exponent also show a 
relatively steep density exponent whereas a steep density expo- 
nent itself does not necessarily require a temperature exponent. 
These internal relations should be explained from the physical 
nature of the clouds. 

There are some pecuUarities concerning four massive cores. 
The fit of the second component of W49A needs a large inner 
region with constant parameters or a very shallow decay of the 
density profile. Here, the main information that we get from the 
line profiles is the central density whereas we can hardly con- 
strain the density exponent in the power-law range. Regarding 
the error bars of the parameters for Serpens we see that the 
the density structure of this relatively nearby cloud is well con- 
strained by its observed size whereas the temperature structure 
is relatively poorly determined. In p Oph it was not possible to 
fit the line profiles simultaneously with the size of the small- 
est core resolved in the observations. This can be explained by 
the clumpy structure of the whole region where the excitation 
of the core at the central position cannot be treated separately 



10 Ossenkopf et al.: Massive core parameters 

Table 7. Direct fit parameters and their likely ranges obtained in the minimisation 
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Table 8. Resulting cloud parameters for S255 using either three 
or four transitions in the fit. 



Fig. 3. Best fit models to the S255 observations fitting only the 
three lower transitions (solid Unes) or all four transitions (dot- 
ted line). 



from the other clumps. Thus we have used for the core model a 
size which is six times larger than the smallest resolved clump 
and contains most of the strong emission observed. To fit the 
profiles and spatial extent of Serpens and SI 06 we have to as- 
sume a relatively large combined abundance (see Sect. |l]) and a 
low hydrogen density in the clumps. This peculiarity could be 
removed when assuming some foreground CS that increases 
the apparent size of the massive core in CS 2-1 but does not 
contribute to the excitation in the core. 



5. Discussion 

5. 1. Restriction by ttie observed transitions 

The main constraint to the core parameters that we can de- 
rive is set by the molecule and the transitions observed. They 
are only sensitive to a relatively narrow density range. We 
can test the limitation of the fits provided by the restriction to 
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three line profiles by comparing the results obtained for cores 
where we have additional CS 10-9 profiles. With the data from 



Plume et al. (1997) for S255 we investigate how much infor- 
mation is lost due to the lack of CS 10-9 data in most cores. 
Fig. ^ shows the resulting best fits to the central line profiles 
in S255 when either all four lines are fitted or only the infor- 
mation from the lower three transitions is used. Tab. || lists the 
corresponding model parameters from the fits. Isothermal mod- 
els provided good fits to the data and the derived physical pa- 
rameters are almost identical except for the inner radius. The 
additional information from the CS 10-9 transition can set a 
smaller limit here corresponding to the higher central densi- 
ties. When predicting the 10-9 line data from the best fit of the 
three other lines the intensity is too low by only about 20% 
(Fig. ||). Hence, we expect reliable results also for those cores 
where only three lines are observed but it would be favourable 
to add information from higher transitions for a better resolu- 
tion of the densest inner region. 

Moreover, a fourth line will reduce the error bars of the pa- 
rameters. In Tab. ^ we compare the error obtained for the S255 
observations using either the three- or the four-lines fit. As dis- 
cussed above the possible range of the inner and outer radii is 
limited only in one direction. Thus, the fourth line mainly re- 
duces the uncertainty of the inner radius. It hardly influences 
the error of the outer radius and the density at the inner radius 
but it also reduces the error of the density exponent, the tem- 
perature and the velocity structure. Thus the inclusion of addi- 
tional lines in the model fits would also give a better constraint 
of the parameters derived. 
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5.2. Comparison of the methods 

Tab. |l^ shows quantities characterising the global properties of 
the clouds computed from the fit parameters in Tab. ^ The av- 
erage density in the second column is given by the cloud mass 
within the outer radius. We see the strong discrepancy between 
the average density (Tab. |l^ and the central clump density 
(Tab. ^ reflecting a very inhomogeneous structure with low 
volume filling factor of dense clumps. 

In the column 3 and 4 of Tab. |l^ one can compare the av- 
erage column density towards the centre in the nonlocal model 
with the column density computed from the escape probabil- 
ity model. Here, we have used the molecular column densities 
from Tab. |] assuming a CS abundance of 4 10~^. We find an 
agreement within a factor of about 1.5, despite the completely 
different analysis applied, except for W49A, W51, and W3 
where the column density from the escape probability model is 
more than a factor of two higher and S235B where it is lower 
S235B, W49A, and W3 are the smallest sources in our sample 
unresolved even in the CS 2-1 beam. Here, the beam filling 
factors used in the escape probability model are uncertain so 
that they may be responsible for the difference. For W51A, no 
simple explanation for the difference is obvious. It is however, 
by far the most massive core in our sample so that it might be 
somewhat pecuhar from that point of view. In general, we find 
that the escape probability analysis provides a reasonable deter- 
mination for the column density when we have a good estimate 
of the beam filling factor It fails to derive correct densities or 
sizes. 



5.3. Core masses 

Columns 6 and 7 in Tab. |lo| shows the core masses computed in 
two different ways. Column 6 gives the integrated mass of the 
model cloud assuming the smallest fitting outer radius and the 
maximum possible central radius. The influence of the inner ra- 
dius on the total mass is negligible, but the uncertainty from the 
lack of information on the outer radius has to be kept in mind. 
Increasing the amount of virtually invisible material around the 
core by increasing the outer cloud radius can easily increase the 
total mass by more than a factor 10. As the mass computation 
relies on the knowledge of the CS molecular abundance (Sect. 
|l]), the resulting values are to be changed if the true abundances 
deviate from the assumed value of 4 10^^. 

Column 7 contains the core virial mass assuming equipar- 
tition of kinetic and gravitational energy in a homogeneous 
spherical cloud. We used the central CS 2-1 line profile and 
the size of the cloud visible in this transition to estimate the ve- 



locity dispersion in the line of sight and the radius. Following 



Lang (1980) we obtain the virial mass by 



Mvir = OMSSAv^DVAaAS 



(3) 



where Av is the FWHM of the line (Tab. ||), D the distance 
of the source, and the Aa and AS the FWHM of the source in 
declination and rectascension in arcmin from Tab. ^. 

For all clouds except p Oph A the agreement between the 
two masses falls within a factor of two. The behaviour of 
p Oph A is due to the radius of the excitation model which 
is larger than the smallest resolved core, as discussed above. 
Hence, it provides a larger mass than the virial estimate which 
uses this visible core size. The general good agreement is quite 
amazing regarding the uncertainty of the outer boundary of the 
models. The clouds seem to be virialised and the CS abundance 
estimate holds approximately for all clouds. 

The agreement of the mass from the SimLine fits with the 
virial mass and independent estimates from the literature indi- 
cates that the cores are well confined and our minimum outer 
radius corresponds to a real, relatively sharp boundary for most 
cores in agreement with the results from continuum observa- 
tions of several cores by Andre et al. (1999). Future investiga- 
tions are, however, necessary to confirm this result because the 
nature of virialisation is still not understood and it is therefore 
not clear how much of the "invisible" low density mass would 
contribute to the virial mass. 

5.4. Comparison with other observations 

To judge the reliability of the parameters derived here, we can 
compare them with core parameters obtained independently 
from observations in other tracers and with other telescopes. In 
general they provide only values for few of the cloud parame- 
ters but they may serve as an independent test of our results. We 
cannot include a complete discussion of the literature concern- 
ing the 15 massive cores considered here. Rather we restrict 
ourselves to a few selected observations showing the general 
power and limitations of the method. 

For NGC2024 we have compared our data with results 
from complementary high-resolution observations in Tab. |ll[ 



Mezgeretal. (1992) combined the results of SEST observa- 



tions at 1.3 mm with IRAM 30m continuum maps at 870 pm 
to identify several clumps in NGC2024 and to deduce their 
physical properties from the continuum fluxes. The given 
values correspond to FIR5 falling at our central position. 



Wiesemeyer et al. (1997) used a spherically symmetric con- 



tinuum transfer model to derive the physical parameters of 
FIR5 from 3 mm continuum observations taken with the 
IRA M Plateau de Bure int erferometer combined with the data 
from Mezgeret al. (1992)| and VLA 1.3 cm observations of 
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Table 10. Clump parameters derived from the SimLine results. The fourth column contains the column density from the escape 
probability model (Tab. |^ for comparison and the last column the distance assumed in the mass derivation 
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Table 11. Comparison of the NGC2024 FIRS parameters with 
results from investigations based on other high resolution ob- 
servations. 
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Wiesemeyer et al. (1997) 

" KOSMA CS and C^*S, IRAM 30m CS as discussed in the text 



Gaume et al. (1992) . Depending on the assumed luminosity 
and dust properties they found a range of parameters fitting the 
observed continuum. The last line in Tab. [ll] represents the re- 
sults from the SimLine fit to our data. In contrast to Tab. |l^ the 
central column density given here is not averaged but computed 
towards a central clump to allow a better comparison with the 
dust observations which are able to resolve this clump. 

We see that one cannot reveal the true radius of the 
core from our low resolution observations. Moreover, the 
CS observations cannot trace the same high densities as the 
dust observations so that they provide only a lower limit. 
It is, however, already close to the central density given by 
Wiesemeyer et al. (1997) The mass and column density de- 



rived from our radiative transfer computations agree quite well 
with the values provided by the high-resolution observations. A 
possible explanation for the difference between the gas kinetic 
temperature and the dust temperatures was provided already by 
Schulz et al. (1991) . They performed NH3 andCS observations 
of NGC2024 and obtained temperatures between 35 and 40 K 
at the position considered. Using a two-component dust model 
they demonstrated that these temperatures are also consistent 



are able to derive realistic values for the core parameters with 
a clumpy radiative transfer model even if we are not able to 
deduce the exact object size as we cannot resolve it. 

W51 consists of three compact molecular cores located 
within about one arcmin. Interferometric observations by 
Young et al. (1998) and the combination of line and contin- 
uum measurements by Rudolph et al. (1990) seem to indicate 
collapse of the component W51e2 with a mass of about 40000 
Mq. The mass determined by the SimLine fit is 54000 Mq. 
The FCRAO and KOSMA observations show no signatures of 
collapse as they are probably blurred by our low spatial res- 
olution. Sievers et al. (1991) obtained temperatures between 
20 K and 57 K and Zang & Ho (1997) derived 40-50 K for an 
inner region of about 0.2 pc and 25-30 K for the outer cloud 
based on NH3 observations. We were able to fit the observa- 
tions with an isothermal cloud at 44 K, however cannot ex- 
clude such a temperature structure. Young et al. (1998) used 
an LTE code assuming spherical or spheroidal symmetry to 
simulate the inner 0.2 pc region of W51e2 fitting the observed 
ammonia data. They obtained density gradients of -1.8 to -2.2 
somewhat steeper than in our fit, indicating that a dense cen- 
tral region might be surrounded by an envelope with a flat- 
ter density gradient. Their central densities between 1.5 and 
22 10^ cm^'^ bracket our value of 5 10^. The assumption of a 
smooth medium by Young et al. (1998) results in a central tem- 
perature estimate below 25 K and a steep temperature gradient, 
i.e. significantly lower temperatures than in our clumpy turbu- 
lent model. We have tested this behaviour by trying to fit the 
data without clumping in our model and also got low temper- 
atures below 20 K but quite bad values. Thus the correct 
treatment of the internal clumping is essential for a reliable 
temperature derivation. 

W3(OH) was studied e.g. by Wilson et al. (1991) using 
VLA obser vations of methanol an d OH and by Tieftrunk et al. 
(1998) and Helmich et al. (1996) with single dish observations 
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of ammonia and HDO, respectively. The KOSMA beam cov- 
ers several maser spots and a luminous mm continuum source 
- probably a class object. Wilson et al. (1991) obtained a 
kinetic temperature of the molecular cores of 20 K, whereas 



Tieftrunk et al. (1998) derived 27 K. Both agreed with our total 
mass estimate of about 2000 Mq. The radius of 1.3 pc deter- 
mined by Tieftrunk is somewhat smaller than the value com- 
puted from our model, but clearly within the error bar. From 



the HDO excitation Helmichet al. (1996) conclude dust tem- 
peratures above 100 K at densities between 10^ and 10^ cm~'^, 
with few embedded clumps at 10^ cm^"^. This clump density 
and the column density of about 2 lO^'^ cm~^ agree approxi- 
mately with the values from our spherical model. The temper- 
ature of 39 K determined from the CS observations falls into 
the range discussed but the difference to the methanol and am- 
monia values asks for an explanation. Around the massive core 
there is probably still an extended envelope of gas at low densi- 
ties insufficient to excite the observed CS transitions (Tieftrunk 
et al. 1998). 

W3 was studied in detail by Tieftrunk et al. (1995, 1997, 
1998) using C^'^S, C^'^O, NH3, and continuum observations. 
Our beam covers the two bright components W3 Main and W3 
West. Moreover, the region contains some ultracompact Hll 
regions related to infrared sources. The molecular line emis- 
sion peaks at a position close to W3 West. The line velocity 
of our CS observations at -42 km/s agrees with the velocity 
of W3 West indicating this component as the main origina- 
tor of the observed CS emission. The combination of single 
dish and VLA observations by Tieftrunk et al. (1997, 1998) 
showed extended gas at a temperature of 25-45 K, a density of 
10^ cm~^, an size of about 1 pc, and a total mass of 1100- 
1400 Mq. This corresponds well to the parameters derived 
from the KOSMA observations. We have traced the emission 
to the somewhat larger radius of 1 .4 pc, but within a radius of 
1 pc we get about the same average density of 0.7 10* cm^'^. 
Our mass estimate of 1 100 Mq and the average temperature of 
41 K also agree. The VLA observations showed several very 
compact clumps with a size of 0.02 pc, densities of 10^ cm^'^, 
and Tkin — 250 K. They are not resolvable from our data, but 
correspond to the clumps in the turbulence description and the 
core parameters derived from the radiative transfer model show 
a similar size and density. In the smooth temperature distribu- 
tion assumed in the radiative transfer model we are not able 
to resolve hot spots with 250 K but found the need for an in- 
creased temperature towards the centre. 

From the group of low mass cores, Castets & Langer (1995) 
analysed CS observations of OMC-2 by means of an LVG anal- 
ysis providing Tkin = 24 K, a density of 9 10^ cm^'^ and a 
CS column density of 5.4 lO^'^ cm^^ in agreement with our 
results. They found akeady indications for substructure with 
clump radii of about 0.022 pc and higher densities in observa- 
tions with higher resolution. Our analysis shows densities of at 
least 1.4 10^ cm^-^ at a scale of 0.019 pc. The virial mass of 71 
Mq computed by Castets & Langer is only somewhat smaller 
than our mass estimate of 95 Mq. Recent 1.3 mm observa- 



data by Lis et. al. (1998) reveal even 30 clumps but lower tem- 
peratures of 17 K supporting our approach of the clumpy cloud 
model. 

The comparison shows that different tracers see different 
parts of a cloud corresponding to different physical conditions. 
Results from other authors based on CS observations agree in 
most cases quite well, whereas the results from other tracers 
may considerably differ The relatively large uncertainty in the 
temperature structure that we cannot resolve within our anal- 
ysis asks for additional observations in higher transitions or 
at better spatial resolution. For nearby clouds like NGC2024, 
OMC-2 or p Oph A we get a good agreement with results from 
high-resolution or even interferometric observations, whereas 
for distant massive cores like W49A and W51A there are sev- 
eral open question, especially regarding the temperature struc- 
ture. 

5.5. The physical nature of massive cores 

All massive cores seem to be approximately virialised inde- 
pendent of their internal structure with respect to the number, 
distribution and luminosity of young stars. Although one could 
expect that violent bipolar outflows observed in some cores will 
drastically change the energy balance in the core, the physics of 
the turbulence in the cores seems to be extremely stable guar- 
anteeing a continuous state of virialisation. 

The relatively sharp outer boundary suggested by the mass 
estimates can be interpreted in terms of collapse models. The 
collapse of an isothermal sphere would result in a self-similar 
density distribution without clear boundary whereas our results 
rather tend towards the scenario of a finite-size Bonnor-Ebert 



condensation ( Bonnor 1956 ). The outer boundary is, however, 
not well determined but only set by the mass constraints be- 
cause the radiative transfer model itself cannot exclude a con- 
tinuation of the density structure to larger radii. 

The density exponent of about -1.6 derived for most cores 



tions by |Chini et al. (1997)| show at least 1 1 embedded conden- 
sations in OMC-2 with masses between 5 and 8 Mq and tem- 
peratures between 20 and 33 K whereas 350 fim continuum 



is consistent with several collapse models (see Sect. 4.3) but 
deviations from the exponent for particular clouds up to values 
around -2 have to be explained. 

On the other hand we have seen that simple collapse mod- 
els are not relevant for the massive cores considered here, 
dumpiness is a main feature of all clouds and smooth micro- 
turbulent models are not able to explain the observed lines. In 
agreement with other high-resolution observations we find typ- 
ical clump sizes of 0.01-0.02 pc at least for the nearby cores. 
In massive distant cores the situation might be more complex 
including a hierarchy of clump sizes resulting in a larger uncer- 
tainty of the temperature profile derived from our model. 

6. Conclusions 

We have shown that the careful analysis of multi-line single 
dish observations with a relatively large beam can provide a 
set of information comparable to single-hne interferometric 
observations. From a careful excitation analysis using a self- 
consistent radiative transfer computation it is possible to de- 
duce some sub-resolution information. We can infer clump 
sizes, masses and densities at scales below a tenth of the beam 
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size. However, interferometric observations are necessary to 
determine the exact core geometry including the location and 
number of clumps within a dense core. 

The spherically symmetric radiative transfer code used here 
is able to take into account radial gradients in all quantities and 
internal dumpiness of the cloud. It enables a reliable deduc- 
tion of the physical parameters from line profiles observed in 
sources with a size close to the spatial resolution limit. The 
method allows to analyse similar observations of objects like 
star-forming cores in distant galaxies unresolvable by all to- 
day's means. For a better resolution of the internal temperature 
structure the approach should be combined with sophisticated 
models on the energy balance including the continuum radia- 
tive transfer in the future. Although the simple escape probabil- 
ity analysis gives a reasonable estimate for the column density, 
it fails regarding the density and temperature structure. 

The line analysis shows two essential points: 
i) The main constraints on the structural quantities which can 
be deduced from the observations are set by the tracer. The 
range of densities and temperatures that one can determine 
from the radiative transfer calculations is restricted by the tran- 
sitions observed. In case of the CS 2-1, 5^, and 7-6 lines, the 
covered densities range from about 2 10^ to 10^ cm^'^. With 
additional information from the CS 10-9 the upper limit can be 
extended by another factor 5. The information from the rarer 
C'^'^S isotope cannot extend the density interval but reduces 
the error bars and provides better estimates for the dumpiness 
of the medium. The high resolution observations discussed in 



Sects. 3.2 and 5.4 show that different tracers provide access to 



different types of information whereas the parameters from our 
CS observations agree well with the CS results there. 
ii) Temperature and dumpiness are related quantities. When 
turbulent clumping in the cloud is neglected, the temperature 
determination will necessarily fail. On the other hand does ac- 
curate information on the clumpy structure of a massive core 
help to constrain the temperature structure. Additional obser- 
vations in higher transitions or complementary estimates of 
the dumpiness will help to reduce the uncertainty of the tem- 
peratures. Thus spatial resolution is still essential. For nearby 
clouds we get a good agreement with results from other high- 
resolution observations, but for distant massive cores the tem- 
perature structure is still an open question. 

All massive cores that we have analysed are characterised 
by turbulent dumpiness with typical clump sizes of 0.01- 
0.02 pc. The clouds are approximately virialised and show den- 
sity gradients around -1.6 but with a scatter between -1.1 and 
-2.2. Large parts of the cores follow a constant temperature but 
we must admit a considerable uncertainty in the most inner and 
outer parts. The correlation between the cloud temperature and 
the turbulent line width indicates that related processes should 
be responsible for heating and turbulent driving. 

Future observations of dense cores should focus on differ- 
ent tracers to gain access to additional information which can- 
not be deduced from a single tracer such as CS. As a drawback, 
the full uncertainty of today's chemical models will enter and 
partially limit the interpretation of the observations. 
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Appendix A: The radiative transfer probiem in 
tlie escape probabiiity 
approximation 

A. 1. The general radiative transfer problem 

The physical parameters of a cloud model and the emerging 
line profiles and intensities are linked by the radiative trans- 
fer problem. It relates the molecular emission and absorption 
coefficient at one point to the radiation field determined by 
the emission and transfer of radiation at other locations in the 
cloud. The quantity entering the balance equations for the level 
populations at a point r is the local radiative energy density u 
within the frequency range for each transition: 



u{r) 
u{r, n) 



dft u{r, n) 



4tt 



1 

- dv I^{r,n)<^^{r,n) . (A.l) 

C J-oo 



Here, w(r, n) is the absorbable radiative energy coming from 
direction n and I^{r,n) is the intensity at a given frequency 
within this direction. It is determined by the radiative transfer 
equation 



nVIy{r,n) = -K.y{r,n)Iy{r,n) + e^(r,n) 



(A.2) 



Assuming complete redistribution the profile for the absorption 
and the emission coefficients, (n, s) and €^{n, s) is given by 
the same local line profile For Maxwellian velocity distri- 
butions it is a Gaussian: 



1 



■ exp 



(A.3) 



Here, v is the velocity of the local volume element written in 
units of the frequency v. The frequently used FWHM of the 
distribution is related to a by FWHM= 2-\/ln2cr. 

For a molecular cloud this results in a huge system of inte- 
gral equations interconnecting the level populations and inten- 
sities at all points within a cloud. 

A.2. The escape probability model 

A simple way to avoid the nonlinear equation system is the es- 
cape probability approximation that is widely applied to inter- 
pret molecular line data. It is based on the assumption that the 
excitation, and thus the absorption and emission coefficients. 
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are constant within those parts of a cloud which are radiatively 
coupled. Then the radiative transfer equation (Eq. A..2) can be 
integrated analytically. We obtain for the integrated radiative 
energy density: 



u{n) 



f3{n)hg{n) + (1 - Pin)) 



(A.4) 



where the subscript I denotes line integrated quantities. The 
term /3(n) is the probability that a photon can escape or pene- 
trate along the hne of sight n from the considered point to the 
boundary of the interaction region or vice versa. 



/oo 
di/ $y(r, n) X exp(-T,,(r, n)) 
-oo 

/r 
dsn K;(r)$^(r,n) 
-oo 



(A.5) 
(A.6) 



where the integration path dsn follows the direction n. 

There are two main concepts to define the interaction re- 
gion and thus to compute the escape probability. The first 
one is the large velocity gradient approximation introduced by 
Sobolev (1957) Here, the interaction region is determined by 
a velocity gradient in the cloud that displaces the line profiles 
along the line of sight so that distant regions are radiatively de- 
coupled. When the resulting interaction region is sufficiently 
small it is justified to assume constant parameters. The escape 
probability then follows 



P{n) = [1 - cxp (-tlvg)]/ ThVG 

Hi 



with 



Tlvg 



{nW){nv) 



(A.7) 
(A.8) 



(c.f. Dssenkopf 1997). The radiative energy density in each di- 
rection is determined by two local quantities only: the source 
function S = ei /ki and the optical depth of the interaction re- 
gion tlvq. The observable brightness temperature at the line 
centre is constant over aU regions with the same velocity gra- 
dient 



2kL 



[S ~ /bg) [1 - exp(-2rLVG)] 



(A.9) 



In molecular clouds, the local velocity gradients are un- 
known, however It is generally assumed that the total observed 
line width, which is composed from turbulent, thermal, and 
systematic contributions, can be used as the measure of the ve- 
locity gradient over the cloud size. This approach was applied 
by Illume et al. (1997^ for the massive cores discussed in the 
text. 

We used another method, the static escape probabil- 
ity model. It does not depend on the velocity structure 
but assumes a special geometry of the interaction region. 
Stutzki & Winnewisser (1985) solved the problem for a homo- 
geneous spherical cloud with constant excitation parameters. 
The resulting escape probability is taken to be constant 



P = di^ X exp(--Ti,) 



(A.IO) 



The surface brightness temperature towards the centre of 
the cloud is given by the same expression as Eq. (A..9) when 
we use the line integrated optical depth at the cloud centre r; 
instead of tlvg- It decays with growing distance from the cloud 
centre. Averaged over the whole cloud, the brightness temper- 
ature at the line centre is given by 

Tb ~ ^ (S - hg) {1 - e{-2Ti)) (A. 11) 



with 



2 

e(x) = ^ (1 - exp(-a;)(l + x)) 



This value would be observed with a beam much larger than 
the cloud. 

When the velocity gradient in the LVG approximation is 
computed from the total line width and the cloud size, it turns 
out that both methods agree when applied to observations with 
a small beam towards the cloud centre. Only f or lar ge-be am ob - 
servations, they differ in the functions in Eq. ( |A.9| ) and ( A. 11 ), 
which are either exp(— 2r) or e(— 2r), but result in similar val- 
ues. 

By set ting up a table of beam temperatures from Eqs. ( |A.9[ ) 
and ( A.. 11 ) and comparing the observed line intensities with the 
tabulated values we can derive three parameters from the obser- 
vations: the kinetic temperature Tkin and the gas density 71H2 
providing mainly the source function, and the column density 
of the considered molecules relative to the line width N/ Av 
providing the photon escape probability. 

Appendix B: SimLine - A one-dimensional 
radiative transfer code 

B. 1. The radiative transfer problem 

SimLine solves the line radiative transfer problem discussed in 
Appx. A. 1 in a spherically symmetric configuration by means 
of a A-iteration. The cod e is similar to the concept described 
by Dickel & Auer (1994) but it contains several extensions and 
achieves a higher accuracy from an adaptive discretisation of 
all independent quantities. 



SimLine integrates the radiative transfer equation (A. 2) for 
a number of rays numerically. In spherical symmetry it is suf- 
ficient to consider the propagation of radiation in one arbitrary 
direction which is taken as z here. The integral is computed 
stepwise from Zi_i to Zi 



Iu{p,Zi) = cxp - / dz K^{p,z) 



dz e,y{p, z) exp 



dz'K^{p, z') 



(B.l) 



where t^, is the optical depth at the cloud centre. 



where p denotes the displacement variable perpendicular to the 
z direction. To minimise the discretisation error the integral 
does not use the source function but only the emission and ab- 
sorption coefficients which are linear in the level populations. 
They are assumed to change linearly between the grid points 
and the exact integration formula for a linear behaviour (which 
is not given here, but can be derived straight forward) is ap- 
plied. This approach provides a reduction of the integration er- 
ror to third order The radial grid is dynamically adjusted to 



16 



Ossenkopf et al.: Massive core parameters 



give a maximum variation of the level populations between two 
neighbouring points below a certain limit. In case of strong ve- 
locity gradients additional points are included on the z-scale 
for a sufficiently dense sampling of changes in the profile func- 
tion. The incident radiation at the outer boundary of the cloud 
is assumed to follow a black body spectrum. 

In spherical symme try th e spatial integration of the radia- 
tive energy density (Eq. A. 1 ) can be reduced to 



u{r) 



2tt 






/ dz' / 


rc J 


— r J — oo 


1 

P 


\/ r'^ — z'^ 



The grid of rays tangential to the radial grid is refined by addi- 
tional rays at intermediate p' values to guarantee a sufficiently 
dense sampling on the z' scale The integration uses a cubic 
spline interpolation. 

With the values of the radiative energy density at each ra- 
dial point and for each transition, the system of balance equa- 
tions can be solved providing new level populations. Here, 
a LU decomposition algorithm with iterative improvement 
(Press et al. 1992) is used. The new level populations are used 
in the next iteration as input for the radiative transfer equation. 
The whole A-iteration scheme is solved using the convergence 
accelerator introduced by Auer (1987)| 

Depending on the physical situation the initial guess is ei- 
ther the optically thin limit, thermalisation or the solution of the 
radiative transfer equation using the LVG approximation (Eq. 

A. 4 ). The stability of the local radiation field is used as con- 
vergence criterion. The number of iterations required for con- 
vergence depends strongly on the optical depth of the model 
cloud. For the examples discussed in this paper only about a 
dozen iterations were necessary but other test cases with com- 
plex molecules like water, non-monotonic velocity gradients, 
and high optical depths require several hundred iterations. 

B. 2. The local turbulence approximation 

The turbulence description uses two additional parameters for 
each spatial point: the width of the velocity distribution a pro- 
viding the local emission profile for optically thin lines and the 
correlation length of the macroturbulent density or velocity dis- 
tribution Icon- 

The width of the velocity distribution cr is composed of a 
turbulent and a thermal contribution 



c 



^turb/ 



(B.3) 



where a Maxwellian distribution of turbulent velocities is as- 
sumed. The relation between the FWHM and the variance of 
the turbulent velocity distribution is given by FWHM(?;tuib ) = 
•^8/3 X ln2 (wtm-ij)). The long range variation of the turbu- 
lence spectrum as described by means of a Kolmogorov or 
Larson exponent is simulated by a radially varying turbulent 
velocity dis persion \/ {vt,„.y,) oc r'^. Exponents 7 between 
about 0.1 ( Poodman etal. 1998| ) and 0.7 (Fuller & Myers 
1992) are observationally justified. 



For the local treatment of coherent units in a turbulent 
medium the considered volume element is subdivided into nu- 
merous clumps with a thermal internal velocity dispersion. 
When each clump is characterized by a Gaussian density dis- 
tribution of molecules with about the same velocity n{r) — 
no X exp(— r^/r^j) the effective absorption coefficient at the 
considered velocity for the whole medium is 



Kcff = "ci X Trr^i 



1 — exp(— r) 



dr 



(B.4) 



where rid is the number density of contributing cells and 
Tc\ — y/TT KTci is their central opacity (Martin et al. 1984). As 
the clumps size 7-ci is the length on which the abundance of 
molecules within the same thermal velocity profile is reduced 
by the factor 1/e, we can compute it from the correlation length 
of the velocity or density structure by = ^con- x crth/c. 

When the turbulent velocity dispersion a is at least three 
times as large as the thermal velocity dispersion ath, we obtain 
an effective absorption coefficient 



gcs 



2 A/ \ '^th 

Trr^i X A{tci) x exp 



{1^ - 1^0? 



(B.5) 



with 



-1 /"OO /'TCXp(— 1)^) 

A{t) = ^ / dv 

Vi" J -00 Jo 



i^^^^^dr' (B.6) 



Here, rtgcs is the total number density of clumps. In case of 
incompressible turbulence, i.e. clumping in velocity space, it 
is equal to the reciprocal cell volume. For small values of the 
clump opacity, A{tc\) is identical to Td and we reproduce the 
microturbulent limit. For large t^i, the function A{tci) saturates 
and we obtain a significant reduction of the effective absorption 
coefficient. In case of density clumps, n^ff (i^) is further reduced 
by the filling factor entering Ug^s ■ In SimLine, this is simulated 
by a corresponding artificial reduction of the molecular abun- 
dance. 



B.3. The central Hll region 

To simulate the effect of a central continuum source in the 
cloud, it is possible to assume an HII region in the cloud centre. 
The Hll region is characterised by two parameters, the electron 
density ric and the kinetic electron temperature Tc- 

The absorption coefficient for electron-ion bremsstrahlung 
in the Rayleigh-Jeans approximation is given by: 



3\/27r (47reomo)3c 



Too 



3/2 



InA 



(B.7) 



where it is assumed that the gas is singly ionised and A is given 
by 



/2/cTe 



V7TOC 



3/2 



47r£omc 



T 



4.96 lOM^ - 



3/2 



Hz 



(B.8) 



for To < 3.6 lO^K. The quantities e and rric denote the electron 
charge and mass, c is the vacuum light velocity, and 7 = 1.781 
(Lang 1980). 
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For a thermal plasma, the emission coefficient follows from 
the Planck function 



(B.9) 



In the radiative transfer computations the frequency depen- 
dence of these continuum coefficients is neglected within the 
molecular line width. Within the H I I region, we substitute the 
molecular coef ficients in Eq. ( p3.l| ) by the quantities from Eqs. 
( ) and ( ) so that we locally switch to a continuum radia- 
tive transfer. 



B.4. Computation of beam temperatures 

When the level populations are known, the beam temperature 
relative to the background is computed by the convolution of 
the emergent intensity with the telescope beam. 

Jo Jo 



mb 



2ki 



(B.IO) 



d0 / dppfn^i,{p,(t)) 



The intensity /^(p) is the value on the cloud surface I^{p) = 
I„{p, V^cioud ~ P^)- We assume a Gaussian profile for the 
telescope beam 



/mb(p, 0) = exp 



-(p-Poffsot)^(l + (/)2) = 



(B.ll) 



The projected beam width is computed from the angular width 
by cr,„b = 7r/648000Dcrnib["] = 2.912 IQ-^ FWEMH 
where D is the distance of the cloud. The program computes a 
radial map with arbitrary spacings. 

B.5. The general code design 

The design of the code is directed towards a high accuracy of 
the computed line profiles. All errors in the different steps of 
the program are explicitly user controlled by setting thresholds. 
All discretisations necessary to treat the problem numerically 
are performed in an adaptive way, i.e. there is no predefined 
grid and all grid parameters will change during the iteration 
procedure. The system of balance equations is truncated when- 
ever the excitation of all higher levels falls below a chosen ac- 
curacy limit. 

Furthermore, the code was pushed towards a high flexibil- 
ity, i.e. the ability to treat a very broad range of physical pa- 
rameters with the same accuracy and without numerical lim- 
itations. The systematic velocities e.g. may range from to 
several times the turbulent velocity and the optical depths may 
vary from negative values for weak masing to values of several 
thousands. 

The program is not optimised towards a high speed. Other 
codes with lower inherentaccuracy may easily run a factor 10 
faster and further improvements are possible. Nevertheless, the 
code is suitable for an interactive work even on a small PC with 
execution times of a few seconds for the models considered in 
this paper. 
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